Posterior maximization
======================

The mode step of estimation -- maximising the posterior kernel (or, with flat
priors, the likelihood) -- is run by an *optimizer*. RISE ships several and lets
you plug in your own; the same machinery is used by maximum-likelihood and
Bayesian estimation, by :doc:`optimal simple rules <../DSGE_capabilities/OSR/OptimizedSimpleRules>`,
and by :doc:`indirect inference <Indirect Inference>` (all of which reduce to
"minimise an objective over a bounded parameter vector").

Choosing an optimizer
---------------------

The optimizer is selected with the ``optimizer`` option, either when building
the model or when calling ``estimate``::

	m    = set(m, 'optimizer', 'fmincon');           % the default
	mest = estimate(m, 'data', db, 'priors', priors, 'optimizer', 'wcsminwel');

To pass options through to the chosen optimizer, give a cell ``{name, opts}``::

	opt          = optimset('fmincon');
	opt.MaxIter  = 2000;
	opt.MaxTime  = 3600;                             % seconds
	opt.TolFun   = 1e-8;  opt.TolX = 1e-8;
	opt.Display  = 'iter';
	mest = estimate(m, 'data', db, 'priors', priors, 'optimizer', {'fmincon', opt});

The common controls (``MaxIter``, ``MaxFunEvals``, ``MaxTime``, ``TolFun``,
``TolX``, ``Display``, and -- for population methods -- ``MaxNodes``) are
honoured across optimizers; ``estimate`` also reports which optimizer was used
in ``print_estimation_results``.

Bundled optimizers
------------------

- MATLAB's ``fmincon`` (the **default**; SQP), and ``fminunc`` / ``fminsearch``
  -- the unconstrained ones are used through bounded wrappers
  (``fminunc_bnd`` / ``fminsearch_bnd``) so that the parameter bounds are
  respected. ``patternsearch`` is available if you have the Global Optimization
  Toolbox.
- ``wcsminwel`` -- a robust quasi-Newton optimizer (Sims's *csminwel*, with a
  restart heuristic) that copes well with awkward likelihood surfaces.
- ``wnewrat`` -- a Newton-type optimizer with a numerically computed Hessian.
- ``wgmhmaxlik`` -- a Metropolis-Hastings-based optimizer that explores the
  surface while climbing it (useful as a robust pre-conditioner before a local
  method).
- ``bee_gate`` -- the artificial bee-colony algorithm, a population-based global
  metaheuristic; see also
  :doc:`../StochasticGlobalOptimization/Stochastic Global Optimization` for the
  broader family of nature-inspired global optimizers.
- a simulated-annealing optimizer.
- ``blockwise_optimization`` -- optimises the parameters block by block, which
  helps when there are many of them.

These live in the ``m/optimizers/`` folder of the toolbox.

Quick reference
~~~~~~~~~~~~~~~

.. list-table::
   :header-rows: 1
   :widths: 18 12 14 12 10 34

   * - Optimizer
     - Family
     - Scope
     - Uses Hessian
     - Parallel
     - Typical use case
   * - ``fmincon``
     - SQP / interior-point
     - local
     - approximated
     - via ``UseParallel``
     - default; smooth, well-scaled likelihoods
   * - ``fminunc_bnd``
     - quasi-Newton (BFGS)
     - local
     - approximated
     - no
     - smooth, no constraints needed beyond bounds
   * - ``fminsearch_bnd``
     - Nelder-Mead simplex
     - local
     - no
     - no
     - very small problems; non-smooth objectives
   * - ``patternsearch``
     - direct search (GPS / MADS)
     - local
     - no
     - via Global Opt Toolbox
     - non-smooth / noisy objectives (requires Global Opt Toolbox)
   * - ``wcsminwel``
     - quasi-Newton (Sims) + restart
     - local
     - returned at end
     - no
     - robust on awkward DSGE likelihood surfaces; common default
       choice after ``fmincon``
   * - ``wnewrat``
     - Newton with numerical Hessian (Ratto)
     - local
     - numerical
     - no
     - mode polishing after a global pass
   * - ``wgmhmaxlik``
     - MH-driven climb (Adjemian)
     - local-with-exploration
     - no
     - no
     - pre-conditioner before a local method on multi-modal
       surfaces
   * - ``bee_gate``
     - artificial bee colony
     - global (population)
     - finite-difference at end
     - yes
     - older global default; superseded by the ``+globalopt`` family
   * - ``rise_lshade`` / ``rise_jso`` / ``rise_cma_es`` / ``rise_de`` /
       ``rise_abc`` / ``rise_aco`` / ``rise_bbo`` / ``rise_agsk`` /
       ``rise_mads``
     - ``+globalopt`` family
     - global (population, or deterministic for ``rise_mads``)
     - finite-difference at end
     - yes (per algorithm)
     - global mode search; see
       :doc:`../StochasticGlobalOptimization/Stochastic Global Optimization`
       for the per-algorithm philosophy and hyperparameters
   * - simulated annealing (``sa.gmhmaxlik``)
     - thermal random walk
     - global
     - no
     - no
     - robust, slow; legacy
   * - ``blockwise_optimization``
     - dispatcher
     - same as inner
     - same as inner
     - same as inner
     - wraps any of the above to optimise parameters block-by-block;
       triggered automatically when ``estim_blocks`` is non-empty

Options accepted across optimizers
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Every bundled optimizer reads the same MATLAB-style fields from its
``options`` struct where they apply:

.. list-table::
   :header-rows: 1
   :widths: 20 55 25

   * - Field
     - Meaning
     - Default
   * - ``MaxIter``
     - iteration cap
     - ``1000``
   * - ``MaxFunEvals``
     - total function-eval cap
     - ``inf``
   * - ``MaxTime``
     - wall-clock cap (seconds)
     - ``inf``
   * - ``TolFun``
     - objective tolerance
     - ``1e-6``
   * - ``TolX``
     - step tolerance
     - ``1e-6``
   * - ``Display``
     - ``'iter'`` / ``'on'`` / ``'off'``
     - ``'off'``
   * - ``MaxNodes``
     - population size (population methods only)
     - algorithm-specific
   * - ``ObjectiveLimit``
     - early-stop floor for the objective
     - ``-inf``

Algorithm-specific hyperparameters (``F``, ``CR`` for ``rise_de``;
``sigma0`` for ``rise_cma_es``; ``limit`` for ``rise_abc``; etc.) are
not part of this common surface. They are documented per algorithm in
:doc:`../StochasticGlobalOptimization/Stochastic Global Optimization`
and can be threaded through ``estimate`` by writing a thin wrapper that
closes over the desired ``hyper`` struct.

Running several optimizers in turn
----------------------------------

A common strategy is to locate the right basin with a global search and then
polish the mode with a fast local method. RISE supports running optimizers one
after another -- each starting from the previous one's solution; see
``help dsge/estimate`` for the exact syntax.

User-defined optimizers
-----------------------

``optimizer`` may also be a function (a function name or a handle) implementing
RISE's optimizer interface -- roughly

.. code-block:: matlab

	[x, f, exitflag, output] = myOptimizer(objfun, x0, lb, ub, options)

where ``objfun(x)`` returns the value to minimise, ``x0`` is the starting
point, ``lb``/``ub`` the bounds, and ``options`` the (RISE-augmented) options
struct. Any of the bundled optimizers (e.g. ``wcsminwel.m``) is a usable
template; see :doc:`../DSGE_capabilities/Extending DSGE`.

Interface specification
~~~~~~~~~~~~~~~~~~~~~~~

The exact signature the estimation engine calls, lifted from
``@dsge/estimate.m``::

    [xfinal, ffinal, exitflag, H] = optimizer(fh, x0, lb, ub, options, varargin)

Inputs:

- ``fh`` -- function handle. The engine calls ``fval = fh(x)`` where
  ``x`` is a column vector inside ``[lb, ub]``. RISE wraps the posterior
  kernel here, so ``fh`` returns a scalar to be **minimised** (the
  *negative* log posterior). The wrapper handles all the bookkeeping
  (parameter assignment, model resolve, filter, prior, penalties for
  general restrictions); the optimizer just sees a scalar objective on
  a box.
- ``x0`` -- starting point as a column vector, dimension ``n``.
- ``lb``, ``ub`` -- column vectors of lower and upper bounds, also
  dimension ``n``. RISE guarantees ``lb <= x0 <= ub``.
- ``options`` -- struct with the MATLAB-style fields listed in the
  *Options accepted across optimizers* table above. Optimizers should
  read what they understand and ignore the rest.
- ``varargin`` -- the trailing entries of the ``{name, ...}`` cell the
  user passed under ``'optimizer'``. For example,
  ``estimate(m, 'optimizer', {'fmincon', 'MaxFunEvals', 1000})``
  forwards ``{'MaxFunEvals', 1000}`` to the optimizer as ``varargin``.

Outputs:

- ``xfinal`` -- column vector of length ``n`` with the best parameter
  vector found. Must satisfy ``lb <= xfinal <= ub``.
- ``ffinal`` -- scalar value ``fh(xfinal)``.
- ``exitflag`` -- integer, in the style of MATLAB's optimizers
  (``1`` = converged successfully, ``0`` = budget exhausted, negative =
  failed).
- ``H`` -- ``n x n`` estimate of the Hessian at ``xfinal``. May be a
  finite-difference Hessian or one produced by the algorithm itself.
  Returning a positive-definite ``H`` is helpful but not required; if
  the optimizer cannot produce one, return ``[]`` and the estimation
  engine falls back to a finite-difference pass.

RISE-added options
~~~~~~~~~~~~~~~~~~

The ``options`` struct is a plain MATLAB-style optimisation options
struct. RISE itself does not inject extra fields into it. What RISE
does add to the *estimation* call is a separate option,
``estim_blocks``: when non-empty, the engine routes the user-chosen
optimizer through ``blockwise_optimization``, which calls the
optimizer once per parameter block. The optimizer sees a *smaller*
``x0`` / ``lb`` / ``ub`` corresponding to one block at a time, but its
own interface is unchanged.

Manual stopping is honoured by every optimizer in the ``+globalopt``
family (``rise_lshade`` etc.): drop a file named
``ManualStoppingFile.txt`` in the working directory and the next
iteration boundary returns the current best.

Minimal worked example
~~~~~~~~~~~~~~~~~~~~~~

A trivial user-defined optimizer that wraps ``fmincon`` and prints a
banner before delegating. The pattern is the same for anything more
elaborate::

    function [xfinal, ffinal, exitflag, H] = myOptimizer(fh, x0, lb, ub, options, varargin)

    fprintf('myOptimizer: %d parameters, MaxIter=%d\n', ...
            numel(x0), options.MaxIter);

    opt = optimoptions('fmincon', ...
        'MaxIterations',           options.MaxIter, ...
        'MaxFunctionEvaluations',  options.MaxFunEvals, ...
        'OptimalityTolerance',     options.TolFun, ...
        'StepTolerance',           options.TolX, ...
        'Display',                 options.Display);

    [xfinal, ffinal, exitflag, ~, ~, ~, H] = ...
        fmincon(fh, x0, [], [], [], [], lb, ub, [], opt);

    end

Plug in by name (the function must be on the MATLAB path) or by
handle::

    m = estimate(m, 'data', db, 'priors', priors, 'optimizer', @myOptimizer);

Any of the bundled wrappers (``wcsminwel.m``, ``wnewrat.m``,
``wgmhmaxlik.m``, ``bee_gate.m``) is a fuller worked example of the
same pattern, including the translation between RISE option names and
the underlying solver's argument list.